A Secular Reiativistic Model For Solar System's Numerical 
Simulations 
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ABSTRACT 

Using Gauss' averaged equations, we compute the secular reiativistic effects gener- 
ated by the Sun on the argument of the perihelion and the mean anomaly of an orbit. 
Then we test different alternative simpler models that have been proposed to repro- 
duce the secular reiativistic effects in the orbital elements. Generally, models introduce 
artificial perturbations that are velocity-independent but that depend on the helio- 
centric distance. If these perturbations are set as an impulse in a constant timestep 
integrator, when the particle approaches perihelion the generated impulse could be 
very strong and badly sampled, originating a spurious orbital evolution. In order to 
overcome this setback, we propose two new models based on a constant, distance- 
independent, perturbation. With these models we obtain the correct secular drift in 
the argument of perihelion and the expected secular orbital evolution is reproduced. 
We also discuss with some detail the secular effect gen erated on the mean anomaly by 
different models. This work is an expanded version of lVenturini fc Gallardd (|2010h . 

Key words: Relativity - methods: numerical - celestial mechanics - comets - aster- 
oids 
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1 INTRODUCTION 

The problem of computing reiativistic effects in planetary systems and in particular in the Solar System can be modeled 
by a perturbation to the Newtonian acceleration and is of increasing importance for future study of low perihelia and low 
semimaj or axis populations. For ou r planetary system, the relevant reiativistic effects are those generated by the Sun (see for 
example Benitez fc Gallardd 200ct ). It is a usual practice in textbooks to show the reiativistic effect in a planet's perihelion 
starting from a simplified problem that gives rise to a unique radial perturbation to the Newtonian potential conserving 



starting from a simplified problem that gives rise to a uni que radial perturbat ion to tne INewtoman potential c onserving 
the angular momentum. But, a more rigor ous analysis jBrumberj 1991 ; Shahid-Saless &: Yeomans 1994 ; Bertotti et al. 20031 ') 

20031 ') gives rise also to a transverse component and the expression for 



intended to agree with the IAU standards ( Soffel et al 
the reiativistic perturbation becomes: 
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where fi = k m© and r, v are heliocentric ( Anderson et alj|l97a : iQuinn et al.l ll991). This expression, that has no normal 



component, is valid when considering the reiativistic effects generated only by a spherically symmetric and non rotating Sun 



When the r otation of the Sun is considered a new effect called gravitomagnetic Lense-Thirring effect appears (|SoffeJll989l : 



Ioridl2005al h Fu rthermore, if all bodi e s have reiativistic contributions a more complex expression should be used as explained 



, for example, Benitez fc Gallardd 1 20081 ). The time in Eq. |T} should be considered as the coordinate time; in order to 



compare wit h observables it is necessary to transform, for ex ample, to Barycentric Coordinate Time (TCB) or Terrestrial 
Time (TT) (|Shahid-Saless fc Yeoman] 1 19941 : ISoffel et al.ll2003r ). 
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In order to avoid either the computational cost or the non sympletic form of Eq. Q, some simpler models that only 
depend on r have been proposed which mimic or reprodu ce the correct relativistic shift in the argument o f the perihelia of 
the bo dies under the gravitation al effect of a central mass jNobili fc Roxburgh|[l98rj ; ISaha fc Tremainjl992h . But, as pointed 
out by Saha fe Tremainel (1994), these models correctly recover the perihelion motion but fail in reproducing the evolution of 
the mean anomaly, M, point that we will discuss in this paper. 

These r-dependent secular models have a disadvantage when used to compute high eccentricity low perihelion orbits in 
constant time-steps algorithms of the type "advance-impulse-advance", as is the usual case in sympletic algorithms. In this 
case, the relativistic perturbation, when computed as an impulse, generates spurious results near the perihelion where the 
relativistic perturbation is stronger and poorly sampled. This is because the expression for the acceleration has a term that 
goes with 1/r 3 and for low r values this impulse can be excessively high. For example, in the Sol ar System the relativisti c 
perturbation for an asteroid at r < 0.1 AU is greater than the gravitational perturbation by Jupiter. ISaha fc Tremaine (1994), 
in a sympletic method for planetary integrations, incorporate part of the relativistic correction ([TJ in the "advance" part of 
the algorithm but an acceleration proportional to 1 /r s remains as part of the "impulse" . 

I n this paper we first review the relativistic effects on orbital elements with special attention to mean anomaly (jRubincam 
19771 ; ICalura et ai1ll997l ; lloriol 120071 ) . Then, we propose a very simple model that correctly reproduces the evolution of the 



argument of the perihelion, even for high eccentricity orbits, and which is very useful for sympletic or, more generally, constant 
timestep algorithms. 



2 TIME AVERAGED GAUSS PLANETARY EQUATIONS 



In order to compute the long term variations in the orbital elements produced by the perturbations it is a standard procedure 
to use the Gauss equations of planetary motion ( Brumber3ll99ll : Murray fc Dermottlll999l : lBeutlerll2005h . As we are interested 
in the secular evolution that these corrections generate, we calculate, in the first place, the Time Averaged Gauss Planetary 
Equations (TAGPE). To do so, we have to insert the radial, transverse, and normal components of the perturbing force (R, 
T and N respectively) into the Gauss equations, evaluate them on an unperturbed keplerian ellipse, integrate the equations 
over one orbital revolution (making use of the change of coordinates dt = \df) and divide them by the orbital period. So the 
mean variation rates are the following: 
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and / is the true anomaly. 



l + e cos / 

These expressions are valid assuming the orbital elements remain constant over an orbital period which is not strictly the 
case. For example, there exist very small amplitude periodic variations in the orbital elements generated by Eq. |T]) ( Richardson 
1988) that have negligible effects on the Gauss equations except for the mean anomaly which is very sensitive to the variations 
of a. At Fig. [T] it is shown the variation of da/dt in one orbital period for Mercury due to Eq.|T}. As the semimajor axis is 
not constant during one orbital period, the equation for < M > should be considered as the secular relativistic effect on a 
body with a semimajor axis which has a mean value given by a. 



3 RELATIVISTIC EFFECTS DUE TO THE CENTRAL STAR 

The relativistic acceleration generated by a non rotating central massive body given by Eq. |T} can be decomposed in a radial 
(R) and transverse (T) components: 
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Figure 1. Daily variations in semimajor axis generated by Eq. JT} as deduced from the Gauss equations for planet Mercury. The 
relativistic perturbation induces an oscillation in the semimajor axis. 
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(see also Damour fc Daruellel 19851 ). Substituting in the TAGPE we can compute the secular variations of orbital elements 
produced by this acceleration: 



(4) 




< M >— 



(5) 



where units are radians per day. 

Eq. @ is the well known relativistic effect on the argument of the perihelion and, for small ecce ntricities, Eq. ([S| coincides, 
for example, with the approximate expressions for the secular drift in M given by llorio ] (I2005bf) and with the secular drift 
in M generated by the relativistic model for low eccentricities used by Vitagliana (1997). The exact formula, valid for all 
eccentricities, is Eq. ([5]). We have checked th is performin g some numerical integ rations with and without relativistic effects 
using the packages Mercury (|Chamberal 19991 ) and Evorb ( Fernandez et alj|2002f ) of a system composed by only the Sun and 
a massless particle with the semimajor axis of Mercury and varying e from to 0.9. (see Fig. |5J. Some caution must be 
taken in computing < M > by means of comparing two numerical integrations. Because of the oscillations on the semimajor 
axis generated by the relativistic correction, a particle integrated with the relativistic model is a particle that during the 
integration, on average, had a mean semimajor axis a 7^ ao being ao the initial value. Therefore, to compute the shift on M 
properly, one must compare the relativistic result with a classical integration in which the particle has a semimajor axis given 
by a, not ao. Otherwise one would be comparing integrations of particles with different semimajor axes, and consequently, 
the obtained AM = M re i — M c i a would be meaningless. 

More problems arise when one wants to take into account the drift on mean anomaly of a specific real body in a 
numerical integration. The uncertainty Aa on the determination of the semimajor axis of real bodies generates an error in the 
mean motion ( An ~ ^/u/a 5 Ao ), which produces an uncertainty of the same order on the drift of mean anomaly. Setting 
An ~< M > and using Eq. ([5]), one finds that if 

Aa> ^(—^=-2) 

then the uncertainty on the osculating elements generates a secular drift on M greater than the relativistic effect itself. That 
means that for properly taking into account relativistic effect on M it is necessary to know with high precision, at least of 
the order of 10 -8 AU for Solar System's bodies, the osculating semimajor axis. That is the case of the terrestrial planets 
from Mercury to M ars which hav e very well defined semimajor axes. On the contrary, the major planets have uncertainties 
significantly bigger (Pitjcva 2008). The relativistic effect in mean anomaly given by Eq. (0) and checked in our numerical tests 
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Figure 2. Secular relativistic effect on mean anomal y in arcsec/c y for a particle with Mercury's semimajor axis. Full line: formula (0; 
dashed line: approximation for low eccentricity orbits <Ioridl2005bh : dots: numerical results comparing two integrations with and without 
the relativistic perurbation given by Eq. ([TJ. 
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Figure 3. Relativistic effect on mean anomaly (arcsec/cy) according to formula l(5} for a massless particle in the gravitational field of 
the Sun. 

is showed in Fig. as a function of the semimajor axis and perihelion distance. For comparison Fig. 0] shows the effect in cj 
given by Eq. Q. 



4 MODELS THAT MIMIC THE SECULAR RELATIVISTIC EFFECTS 
4.1 Proposed Models 

If we are interested in computing the relativistic effects due to the Sun on a small body, we could just introduce Eq. flj 
into an integrator. The problem is that this acceleration depends on both vectors position an velocity of the particle, so 
the speed of the integrator may be slowed down in order to calculate accurately the vectorial products at small heliocentric 
distances. Moreover, it i s not a simple task to introduce this perturbation in a sympletic integrator, though it can be done 
i Saha &: Tremainel 1994 V In order to overcome this difficulty, some alternative simpler models have been created in the last 
two decades. The relativistic precession of the argument of perihelion is correctly reproduced defining a radial (T = 0) 
acceleration: 

(6) 



( Nobili & Roxbure 
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for mean anomaly. 


Saha & Tremainel 



illl98rj|. Inserting t his R in the TAGPE the exact secular drifts generated by Eq. {1} are recovered except 
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Saha fc Tremainel (|l992h added one more term into ([6j in order to account for both u) and M drifts: 
4 



(7) 
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Figure 4. Relativistic effect on the argument of the perihelion (arcsec/cy) according to formula $4$ for a massless particle in the 
gravitational field of the Sun. 



Now, ins erting this perturbat i on in the TAGPE, Eq. Q and Eq. ([5} are recovered. This apparen t improvement to the original 
model of iNobili fc Roxburgh (1986) is not so at all as pointed out by ISaha fc Tremainel l| 19941 ). The point is that given an 
initial osculating ao the mean a generated by the evolution under Eq. (JTJ) is different from the mean clr generated by the 
evolution under Eq. (0. Then, model (JTJ) will not reproduce the secular evolution of an object with the correct a but with a 
different mean semimajor axis given by clr. In consequence, the model introduces a secular drift in M and no evident progress 
is done in comparison with Eq. JB). This problem can be solved taking appropriate initial conditions as we explain in sec. 4.3. 



4.2 Our r-independent Models 

The corrections exposed above, though computationally better than Eq. {TJ because the v-dependance is eliminated; still have 
the problem that near perihelion the perturbation can be high enough to introduce numerical errors in constant time-step 
integrators. 

One way of overcoming this difficulty is to look for a constant r-independent relativistic correction. To accomplish this 
task, we went back to Eq. (O, and assuming that the short period terms in the orbital elements do not affect the secular 
evolution of them, we time averaged the relativistic perturbation produced by a central body, obtaining: 

< Ar >= ;£= u e 

c 2 aV(l-e 2 ) 3 

where u e is the versor pointing to the pericenter. The fact that this averaged acceleration is a constant along u e , gave us 
the idea of trying with a constant perturbation in this direction. In order to obtain the same drift in ui as Eq. @ using the 
TAGPE it is necessary to introduce a factor 2 and in terms of R and T can be written as 



R = cos r 

c 2 aV(l-e 2 ) 3 

T = sin / (8) 

c 2 aV(l-e 2 ) 3 

Inserting these expressions into the TAGPE, the elements a, e, i, have zero variation and the mean variation in the u 
is recovered, but for M we obtain: 



<M>=-Mi+£!) 



a 5 (l-e 2 ) 3 



which is different from the real effect given by ((5}. Even so, one can also use a simpler constant radial perturbation (T — 0) 
that generates the expected variation in the argument of perihelion: 

R= 3 JL= (9) 

c 2 aV(l-e 2 ) 3 

Again, from the TAGPE, the elements a, e, i, f2 have zero variation and the mean variation in the w is recovered but for mean 
anomaly we obtain a different expression: 
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Figure 5. Results for planet Mercury from a Solar System's relativistic numerical integration. Evolution of eccentricity and inclination 
according to models JT}, 10 and lUpl. The numer i cal res ults of the all three different models coincide. A comparison with a non relativistic 
integration can be found in Be nitez fe Gallardo! ( 2008 ) . 



<m>=--LJ / 3 (io) 

c 2 y a 5 (l - e 2 ) 3 

but which is coincident with ((5| for e — > 0. Therefore, with both constant models the drift in cj is recovered but there is 
a discrepancy for the drift in M. This is irrelevant in numerical simulations because, as we have explained, the predictions 
for M given by models cannot coincide with the true relativistic effect if the same set of initial conditions for all models are 
taken. 

The last model that we propose (the constant radial acceleration) is somehow better than the one along u e because is 
computationally less demanding and also, in the case of low eccentricities, both Eq. ([5]) and Eq. (|10l) give the same result. 
Moreover, since model ((9]) only depends on (a, e) it is not necessary to evaluate the magnitude of the perturbation in each 
time step. The most important point is that model ((9| maintains the precision of the numerical integration even for very 
small perihelion orbits while the original Eq. fTJ needs a strong reduction in the step size. 



4.3 On Initial Conditions 



Initial conditions should be used according to the theory they were determined. For example, the model of I S aha fe Tremaine 



(1992) can be used confidently for generating precise ephemeris only if the initial conditions were obtained adjusting obser- 



vations to this mo del. That is the u nderlying idea of the very accurate method for computing ephemeris of low eccentricity 
orbits proposed by Vitaglianol (1997). But, it is not possible to follow the exact evolution of the mean anomaly if there is no 



consistence between the initial conditions (which they also should be very precise) and the model used. Then, any of the mod- 
els we have analyzed are valid to follow the secular evolution of an orbit and, in particular, our constant radial perturbation 
model is computationally more convenient. 

For illustrative purposes, in order to compare the variations of the orbital elements produced by the different models given 
by Eq. |l|), Eq. (0 and Eq. (J9j> , we integrated numerically our planetary system using Evorb for one million years. Results 
for planet Mercury are shown in Fig. [5] and Fig. [6] These figures corroborate what we have already stated by means of the 
Gauss equations: the secular relativistic evolution of all orbital elements is very well reproduced except for mean anomaly (a 



comparison with the classic evolution can be found in lBenitez fc Gallarddl2008r ). And even with strong differences in M (Fig 



[6]bottom), the orbital evolution of the different models is almost undistinguishable. For example, at the end of the integration, 
differences of the order of 10~ 6 AU in a, 1(T 5 in e, 1(T 3 de grees on lj, 10" 4 de grees on Q, and even lower differences for i 
were obtained. Note also that the numerical deviations on M with respect to model flj in Fig. [6] are mostly a consequence of 
using the same initial conditions for integrations that deal with different models. In the case of model ([7]), the drift on M will 
disappear if initial conditions consistent with this model are taken. In the case of model ([9]) some drift will persist. 



5 CONCLUSIONS 

By means of numerical integrations we confirmed the predicted secular drift in mean anomaly given by formula ([5]) , which is 
valid for all eccentricities and was obtained introducing the relativistic effect due to the Sun (Eq. [T]) into the time averaged 
Gauss' planetary equations. Several models are able to reproduce the drift in u but for high eccentricity orbits only the model 
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Figure 6. Top: evolution of Mercury's argument of perihelion according to models JTJ, Q and J9j. Al three models are coincident. 
Bottom: differences on Mercury's mean anomaly between model (O and Jl} (full line) and between model JTJl and {I} (dotted line). In 
spite of these differences the secular evolution is the same. 



bv lSaha fe Tremaind (|1992l ) can follow the exact secular drift in M and only if appropriate initial conditions are taken. For 
secular evolution studies of numerous populations by means of numerical integrations, if no precise ephemeris are required, our 
constant radial independent model is very convenient, specially if low perihelion orbits are involved. For accurate ephemeris 
computations of real bodies, the original formula {1} or even the full relativistic N-body one should be used with a control of 
the integrator's precision in the case of low perihelion orbits. 
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